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ABSTRACT 

We study gravitational lensing with a multiple lens plane approach, proposing a simple analytical 
model for the probability distribution function (PDF) of the dark matter convergence, k, for the 
different lens planes in a given cosmology as a function of redshift and smoothing angle, 6. The model 
is fixed solely by the variance of k, which in turn is fixed by the amplitude of the power spectrum, 
(jg. We test the PDF against a high resolution Tree-Particle-Mesh simulation and find that it is far 
superior to the Gaussian or the lognormal, especially for small values of 9 <C V and at large values of 
n relevant to strong lensing. With this model, we predict the probabilities of strong lensing by a single 
plane or by multiple planes. We find that for 9 ~ 10", a single plane accounts for almost all (~ 98%) 
of the strong lensing cases for source redshift unity. However, for a more typical source redshift of 4, 
about 12% of the strong lensing cases will result from the contribution of a secondary clump of matter 
along the line of sight, introducing a systematic error in the determination of the surface density 
of clusters, typically overestimating it by about 2 — 5%. We also find that matter inhomogenieties 
introduce a dispersion in the value of the angular diameter distance about its cosmological mean. 
The probable error relative to the mean increases with redshift to a value of about 8% for z ~ 6 and 
6 ~ 10 ". 

Subject headings: cosmology: theory — gravitational lensing large-scale structure of the universe 
— methods: analytical 


1. INTRODUCTION 

Gravitational lensing provides a unique tool for study¬ 
ing the cosmological distribution of matter, because it di¬ 
rectly probes the gravitational potential and is indifferent 
to the nature or the physical state of matter. It is, there¬ 
fore, free from the assumptions that can plague other 
dynamical methods for studying gravitational fields. 

Tidal deflections of light rays by structures along its 
path can weakly distort and magnify distant sources. 
These effects, although small on the individual level, 
can be statistically significant for an ensemble of sources. 
Thus, statistical properties of the intervening matter dis¬ 
tribution can be inferred from such weak lensing studies 
(see, for e.g. Blandford & Narayan 1992; Bartelmann & 
Schneider 2001; Refregier 2003, for reviews). This field 
is extremely promising at the present time but is fraught 
with observational challenges, because these subtle ef¬ 
fects have to be very delicately isolated from comparable 
systematic errors. 

Occasionally, the distortion due to the intervening mat¬ 
ter may be strong enough to produce multiple images of 
the source. Strong lensing is extremely sensitive to the 
cosmology and hence, can be used to constrain cosmo¬ 
logical parameters. Small changes in the cosmological 
parameters can significantly alter the nonlinear struc¬ 
ture formation history of the universe and thereby result 
in huge differences in the abundance of compact objects 
which behave as potential strong lenses. This makes the 
probability of any specific type of strong lensing event 
a useful observable for constraining cosmology (see, for 
example, Turner et al. 1984; Narayan & Blandford 1991; 
Kaiser 1992; Cen et al. 1994; Nakamura & Suto 1997; 
Cohn et al. 2001; Keeton 2001; Premadi et al. 2001; Oguri 
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et al. 2002; Li & Ostriker 2002, 2003). For example, the 
statistics of lensed quasars or galaxies lensed by clus¬ 
ters into giant arcs can be used to study the population 
of lenses, their distribution in redshift space, mass, and 
hence the cosmological model. Strong lensing systems 
may also be studied individually, in order to obtain a de¬ 
tailed understanding of the structure and the distribution 
of matter in the lensing object and the source (see, for 
example, Claeskens et al. 2000; Broadhurst et al. 2000; 
Shapiro & Iliev 2000; Cohn et al. 2001; Belokurov et al. 
2001; Patnaik & Narasimha 2001; Inada et al. 2003; Fass- 
nacht et al. 2005; Sluse et al. 2005; Wayth et al. 2005; 
Claeskens et al. 2006). In this paper, we shall be mostly 
interested in the former problem, namely the statistics 
of strong lensing events. 

In the absence of satisfactory analytic tools, the most 
widely used method for studying strong lensing statistics 
has been the analysis of large numerical simulations. The 
method of tracing the path of light rays through a simu¬ 
lated universe evolving in accordance with some specific 
cosmological model is being extensively used to study 
the multitude of lensing phenomena. For example, giant 
arc statistics using ray tracing simulations of clusters in 
a ACDM cosmology are being thoroughly investigated 
(examples include Bartelmann et al. 1998; Wambsganss 
et al. 2004a; Dalai et al. 2004; Ho & White 2004; Li et 
al. 2005; Hennawi et al. 2005a,b, etc.). However, these 
simulations are costly in that they require extensive com¬ 
puting power and time. This also restricts the range of 
cosmological models that can be investigated; in most 
cases only the currently favored cosmological model is 
used. One must not, however overlook the array of im¬ 
portant analytical works from the so-called halo model 
approach (see, for e.g. Li & Ostriker 2003; Oguri et al. 
2003; Oguri & Keeton 2004), which allow accurate com- 
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putation of lensing probabilities for a single lens. 

The goal of this paper is to provide an analytic method 
to address some of the questions pertaining to gravita¬ 
tional lensing statistics that, at the cost of having much 
more limited information available than from detailed 
numerical studies, would have the benefit of allowing any 
definite cosmological model to be investigated. We test 
our results against numerical simulations performed by 
the Tree-Particle-Mesh (TPM) code (Bode & Ostriker 
2003) for a ACDM universe. The method lends itself 
easily to the treatment of lensing under other cosmolo¬ 
gies. 

We consider the fluctuations in the matter density 
along a light cone of a given opening angle. We pro¬ 
pose an analytic method to compute the probability 
that a light beam will encounter a given surface density 
while passing through a slice of matter centered around 
a given redshift. The model is completely defined by 
the variance of the projected density of the matter in the 
slice, smoothed on the angular scale corresponding to the 
beam. We show how the variance of surface mass density 
on the lens can be predicted from the power spectrum in 
Section 2 and develop the analytic model in Section 3. 

Armed with this result, we adopt a multiple lens plane 
approach to study strong lensing. We divide up the line 
of sight distance to the source into contiguous slices of 
equal comoving depth, each of which is chosen to be 
larger than the correlation length so that each slice can 
be considered statistically independent of the others. We 
imagine the matter inside each of these slices to be pro¬ 
jected onto a plane at its center. The probability den¬ 
sity function (PDF) for the surface density on each plane 
is then constructed with our analytical method. The 
matter encountered by the light ray while propagating 
through the universe is then simulated by randomly pick¬ 
ing a surface mass density from the PDF on each slice 
it passes through. In fact, the surface density scaled by 
the so called critical density required for multiple image 
formation, better known as the convergence, is computed 
on each of these planes. We identify different strong lens¬ 
ing scenarios according as whether the individual or the 
additive value of convergences so obtained is in excess 
of unity. We shoot a statistically significant number of 
rays to obtain the probabilities of different events. In 
particular, we study the probability of strong lensing as 
a function of source redshift and angular scale. 

Another important question is the role of auxiliary 
planes aiding the main lens to produce multiple images. 
It is known that some multiply imaged quasars (e.g. the 
quadrupole system B1422+231, see Kormann et al. 1994) 
cannot be well modeled under the thin lens approxima¬ 
tion, which assumes all mass responsible for lensing to be 
concentrated on a mathematical surface of zero thickness 
(Schneider et al. 1992). Also, observations suggest cases 
where galaxies at disparate redshifts appear to have con¬ 
tributed to the lensing event (Tonry 1998; Augusto et 
al. 2001; Chae et al. 2001; Johnston et al. 2003; Fass- 
nacht et al. 2005). In order to study the probabilities 
of such events, we identified, within the strong lensing 
cases, the ones that were caused by a single plane or by 
the concordance of two or more planes and studied their 
probabilities separately. Incidentally, a similar study was 
done in Wambsganss et al. (2004b), who studied the phe¬ 
nomenon at the pixel level (~ \.b") and reported a signif¬ 


icant contribution of auxiliary planes in the strong lens¬ 
ing probabilities. With our method, we can conveniently 
include an element of angular smoothing depending on 
the strong lensing situation and study the probabilities 
as a function of smoothing radius. We treat this problem 
in Section 4. 

We also note that our model works very well even in the 
weak lensing regime (angles > 1'). The statistics of weak 
lensing convergence can be studied from the magnifica¬ 
tion effects of clustered matter which leads to variations 
in the number density and image sizes of galaxies across 
the sky (see, e.g., Jain 2002). Also, shear maps generated 
from weak lensing surveys yield statistical properties of 
the convergence. Thus a good analytic model for the 
convergence can be useful in constraining the properties 
of large scale structure (Barber et al. 2004), although 
the practicalities involved may seem rather daunting at 
the moment (Bartelmann & Schneider 2001). The situa¬ 
tion will improve with forthcoming surveys having high 
signal-to-noise ratios, enabling us to study the small scale 
clustering of matter. 

In Section 5, we study the effect of the chance accumu¬ 
lation of secondary matter along the line-of-sight on the 
estimation of cluster masses and finally, in Section 6, we 
examine the dispersion in the angular diameter distances 
due to inhomogeneous distribution of matter. 

2. VARIANCE IN PROJECTED OVERDENSITIES 

Since the fluctuations in matter density along the line 
of sight parallel the growth of structure in the universe, 
with higher redshifts corresponding to an earlier era in 
the structure formation history, these variances are in¬ 
tricately related, and can, in fact, be predicted from the 
evolution of the matter power-spectrunr. For large scales 
or early times the growth of structure is linear and the 
power-spectrunr is linearly related to its primordial form 
through a transfer function (see Bardeen et al. 1986). 
However, on small scales and late times, the growth of 
structure becomes nonlinear and small scale modes cou¬ 
ple and grow more rapidly than in the linear regime. 
The evolution of the power-spectrunr in this era cannot 
be calculated analytically with sufficient accuracy and 
one has to resort to numerical simulations to find some 
analytic fitting function for the simulated power spec¬ 
trum. This was first done by Peacock & Dodds (1996) 
who used N-body simulations based on hierarchical clus¬ 
tering of matter to produce an analytic fit. Later, Smith 
et al. (2003) proposed a halo model (Seljak 2000; Ma & 
Fry 2000) inspired fitting function which was shown to be 
in excellent agreement with numerical simulations (Bar¬ 
ber & Taylor 2003) in the context of the lensing conver¬ 
gence power-spectrum. We have adopted the Smith et al. 
(2003) fitting functions in this work. In order to calcu¬ 
late the variances in the 2D matter distribution obtained 
by projecting the matter along the line of sight in a red¬ 
shift slice, the 3D power-spectrum has to be related to its 
2D counterpart. This is done by the well known Limber 
approximation (Kaiser 1998). Using this method, we 
show that the 2D variances and consequently, the prob¬ 
ability distribution for surface density can be predicted 
with sufficient accuracy given a cosmological model. 

In the present approach, we divide the matter along 
the line of sight from the observing redshift of zero to a 
redshift of 6.34 into 38 slices, each 160 h _1 Mpc thick, 
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and project the matter in each of these slices onto a 
plane at its center. Therefore we have lens planes at 
80 + j x 160 h -1 Mpc for j = 0, ..,37. Let E(£, z) de¬ 
note the projected surface mass density at a position 
£ = {£i,£ 2 } on a plane at redshift z. We denote the 
same quantity, when smoothed by a Gaussian window of 
angular radius 8 q, by Eg 0 (£, 2 ). We are interested in find¬ 
ing the variance = (<5 2 (£>0o)) in the projected 

overdensities <52 (£) given by, 


5 2 (£; do) = 


Ee 0 (£,z) - (Se 0 


Here and henceforth, the subscript 2 on any symbol is 
meant to assert that it is a quantity depending on 2- 
dimensional variables. In order to obtain 02 ( 60 , z), we 
first use the Limber Equation (Kaiser 1998) to obtain the 
2-D power spectrum for the distribution of E on a plane 
from the nonlinear 3D power spectrum i^ L (k). Let us de¬ 
note the comoving distance to the center of the bounding 
planes of the ith slice by Xi-i and x»> so that the first 
slice is bounded planes at by xo and xi and so on. The 
coordinate of a point inside the ?'th slice can be denoted 
by, 

x(x, 0) = x(z)( 6 i, 02,1), (Xi -1 < X < Xi), 

where 6 1 and 62 are two dimensionless numbers (angles) 
specifying the position of the point on the plane orthogo¬ 
nal to the line of sight. Since we are interested only in the 
projection of the matter inside this slice, the overdensity 
seen at the point (# 1 , 62 ) on the plane of projection is, 

6 2 ( 0 ) = J d X W(x) 6 [x(x, 0 )], (2) 

where IT is a selection function, which, in our case, is 
given by a step function, 

W(x) = I fe-L-x) if (Xi-1 < * < Xi) (3) 

W \ 0 otherwise. y J 

The Fourier space conjugate to 9 is denoted by Z, and 
defined through, 


S 2 (l) = / d 2 0 e a e S 2 ( 0 ). 


The 2D power-spectrum is given by the standard Limber 
approximation (Kaiser 1998; Dodelson 2003), 


r w 2 ( y) 

P 2 (l)= d X ^-^-P NL (l/x,z) 
J X. 


= _I_ r dx PM) (6) 

(y-y ,j 2 / X y 2 ’ W 

A4 —1/ J X =Xi-l A. 

using our selection function. 

By introducing k = l/x, we can rewrite this as, 

1 rWxi-i 

P 2 (l,z) = - -Tj / d kP NL (k,z). (7) 

HX* - Xi-iY Ji/ Xi 

Now, the variance of the projected overdensity, Gaussian 
smoothed at an angle 6 q can be calculated by using, 

°2 ( 0°> z i) = I A2 ( W °W> d X (8) 


where A(l,9o) is the Fourier transform of the 9 space 
Gaussian smoothing window of width $o, and is given 

by, 

A(Z,0 o ) =e~ l 2 e Y 2 . 

With this we have, 


a 2 (e 0 ,Zi) = j^^ J d 2 l e l2e °P 2 {l) 


27r(Xi - Xi- 1) 2 


2 2 r l /xi -1 

d le~ l0 ° dfcP NL (fc,2,;) 

Jl/Xi 


(Xi - Xi- 1) 2 


d l e~ l202 ° 


r^NL, ( 10 ) 


where A^j L = -^pPni J {k)k d ‘ is the dimensionless power 
spectrum. 

The integral above requires the value the non-linear 
power spectrum A§j L (fc) as a function of scale k and red- 
shift 2 . Incidentally, this can be generated using the 
routine “halofit” by Smith and Peacock (see Smith et al. 
2003). 

We test our predictions against a high resolution 
cosmological simulation using the Tree-Particle-Mesli 
(TPM) code (Bode & Ostriker 2003). The simulations 
were performed in a box with a comoving side length of 
L = 320/i _1 Mpc. We used N = 1024 3 = 1,073, 741,824 
particles, so the individual particle mass is to p = 2.54 x 
10 9 ft. _1 Mq. The cubic spline softening length was set 
to e = 3.2 h ^ 1 kpc, producing a ratio of box size to soft¬ 
ening length of L/e = 10 5 . The output was stored at 
19 redshift values out to z ~ 6.4, such that the cen¬ 
ters of the saved boxes matched comoving distances of 
(160+nx320)h _1 Mpc, where n = 0,..., 18. Each box was 
then divided into 9x9 separate square cylinders along its 
length. This was done for three orthogonal projections, 
leading to 243 sub-volumes. Two lens planes were pro¬ 
duced for each sub-volume by bisecting along the line-of- 
sight and projecting the mass in each 160/i -1 Mpc-long 
volume onto a plane. Thus, we have planes at comoving 
distances of (80 + j x 160) h _1 Mpc, where j = 0, ...,37. 
At the highest redshift, each plane has a side length of 
35.6/i -1 Mpc and contains 800 2 pixels, making the pixel 
size ddAh ^ 1 kpc comoving. At lower redshifts, an open¬ 
ing angle of about 20" was maintained by keeping the 
number of pixels constant but progressively decreasing 
the size of the planes. Thus in the lowest redshift box 
the lensing planes are 1.9/i _1 Mpc on a side and the pixel 
size is 2.3/i _1 kpc, which is slightly lower than our spatial 
resolution. Because at lower redshifts the plane does not 
cover all of the corresponding sub-volume, there is a ran¬ 
dom offset perpendicular to the line of sight within each 
sub-volume. In all lensing applications discussed below, 
we will not consider structures smaller in size than 5 pix¬ 
els or three times our spatial resolution. 

In order to test the results against the WMAP normal¬ 
ized (Spergel et al. 2003) ACDM TPM simulations, we 
used the same values of cosmological parameters as used 
by the TPM code, to generate the power spectrum by 
halofit. These parameters are, fl m = 0.3, = 0.7, 

fib = 0.04, dimensionless Hubble constant, h = 0.7, 
<78 = 0.95, and the derived power-spectrum shape param¬ 
eter (Sugiyama 1995) F = fl m hexp —fib ^1 + • 
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Fig. 1.— Top : Comparison of the analytic result for the variance 
a! (z, 6) as a function of redshift z for three different smoothing 
angles 6 , with the data obtained from the quoted TPM numerical 
simulations (Bode Sz Ostriker 2003). The points are the variances 
from the N-body simulations. The lines represent the analytic 
result obtained by using the Limber equation and the nonlinear 
power spectrum of Smith et al. (2003). Bottom : The percentage 
errors in the predicted values of the variance 


We used redshifts corresponding to the center of each 
of the 160 h -1 Mpc thick slices to generate the power 
spectrum for that slice. Then we used the Limber Ap¬ 
proximation and the Gaussian filtering by the desired 
angle [equation (10)], to obtain the variances of the 2- 
D overdensities on the planes at each of these redshifts. 
The results are displayed in Fig. 1. This test shows that 
we are able to predict the variance within each of these 
slices from theory to roughly 1 — 10% accuracy. 

3. THE MODEL 

The probability distribution function for the lensing 
convergence for large angular scales, > 10', can be 
straightforwardly obtained from analytic perturbative 
calculations (Stebbins 1996; Villumsen 1996; Bernardeau 
et al. 1997; Jain & Seljak 1997; Kaiser 1998; Schnei¬ 
der et al. 1998; van Waerbeke et al. 1999). However, 
on angular scales less that 5', the perturbative approach 
fails and non-linear gravitational clustering comes into 
play. Hence, a number of models have been developed 
to move from a linear to the fully nonlinear power spec¬ 
trum, and from there to the PDF of the convergence 
field. On intermediate scales (1 — 10') such analytical 
models have been shown to have good predictive power 
when compared with numerical simulations (see Barber 
et al. 2004, and references therein.). Recently, Kayo et 
al. (2001) analyzed the one- and two-point statistics of 
the three dimensional dark matter distribution using N- 
body simulations and found them to be well described 
by the lognormal model (Kayo et al. 2001). Following 
this, Taruya et al. (2002) have shown that the lognor¬ 
mal model, with certain parameterizations, can fairly 


describe the convergence field for small angular scales 
2' < 0 < 4' and up to a level threshold v ~ 10, where 

v = k/ (k 2 ) . However, satisfactory analytic methods 

to generate the convergence PDF on the ultra non-linear 
scales, < 1', are as yet unknown. Here, we put forward a 
phenomenological model for the PDF of the surface mass 
density and hence, that of the convergence on the lens 
sheets, which seems to work down to few arcseconds and 
is designed to better approximate the high density tail of 
the distribution, which reduces to a distribution akin to 
the lognormal and then to the Gaussian for large angles 
and high redshifts. The model is completely defined by 
the variance of the surface mass density on each sheet 
and thus encodes all higher order statistical information 
in the two point statistics of E. 

We computed the probability density function (PDF) 
of the quantity x = = 1 + 62 on each plane. It is 

found that the PDF is well described by the following 
form, inspired by the lognormal distribution, 


,, , N 

f{x) = — exp 
x 


(In a; + w 2 /2) 2 (l + A/x) 
2^2 


( 11 ) 


where the values of the three parameters, A, u> and N 
are fixed by the following three constraints, 

Normalization : f(x)dx = 1, (12) 

By defn. (x) = 1 : / 0 °° xf(x)dx = 1, (13) 

Observed variance: J 0 ° °(x - (x)) 2 fdx =a\. (14) 

Therefore the model is fixed by a single input, the vari¬ 
ance cr 2 (0, z) of the projected overdensity of at a given 
redshift and a given angular scale of smoothing. Since 
this quantity can be theoretically predicted, as discussed 
in the previous section, the model is essentially fixed by 
the cosmology. 

A search algorithm is developed which, given a value of 
cr 2 , returns the values of the three parameters satisfying 
all above constraints. 

The result is displayed in Fig. 2, where the parameters 
A, lu 2 and N are plotted as functions of <J 2 (x). Note that 
all these parameters are approximately power-laws of cr 2 
for small values of the latter. 

The theoretical model is shown against the PDF ob¬ 
tained from simulations for angular smoothing radii of 
7.5", 52.5" and 1.5' in Figs. 3, 4 and 5, respectively, 
with the lognormal model plotted over each. A detailed 
comparison of the model and other standard PDFs are 
displayed in Fig. 6. 

We would like to stress here that although our model 
has three variables, A, to and N, it is not a three parame¬ 
ter fit in the sense that the three variables are not varied 
to obtain some kind of “best” fit with the data points. 
We have simply put forward a non-standard probability 
distribution which is on the same footing as the Gaussian 
or the lognormal as far as the number of free parameters 
is concerned. All the three distributions have their vari¬ 
ables determined by the fact that they are normalized, 
have unit mean by construction and have a certain vari¬ 
ance predicted by the cosmology. Thus, each of the three 
distributions is essentially fixed by the input cosmology 
and the knowledge of the redshift and the smoothing an¬ 
gle. In Figs. 3, 4 and 5 we have simply overlaid the 
theoretical model on the data points; there is no fitting 
with the data involved anywhere. 
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Fig. 2.— The parameters A, lo 2 and N plotted as a function of 
( t 2 (x ) (cf. equation (11)) 



Fig. 3.— The theoretical model (red line) shown against the 
probability densities (blue circles) obtained from the quoted TPM 
simulation (Bode & Ostriker 2003) for different redshifts. The 
dashed dark blue line is the lognormal model. The smoothing 
angle for this set of data was 7.5 r/ . The color coded arrows show 
the position of the critical surface density E cr it in units of (E), for 
four different source redshifts: 7.0 (black), 5.0 (green), 3.0 (blue) 
and 1.0 (red). 


It is immediately apparent from these figures that our 
model does better than the lognormal or the Gaussian in 
describing the PDF of the surface density, in particular 
it better describes the all important tail of the distribu¬ 
tion which characterizes the low probability but observa- 
tionally important high k events. The lognormal hugely 
underestimates the probability density for just those rare 
events in the tail of the distribution which are important 
for strong lensing. This effect is more pronounced for 
smaller smoothing angles, as evident from Fig. 3, where 
we show, by arrows, the critical surface densities required 
to produce strong lensing, for four representative source 
redshifts. At these values, the probability densities pre¬ 
dicted by the lognormal model are orders of magnitude 
lower than those obtained from the TPM simulation. For 
all smoothing scales, it appears that on the right tail of 
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Fig. 4.— The theoretical model (red line) shown against 
the probability densities (blue circles) obtained from the quoted 
TPM simulation (Bode Sz Ostriker 2003) for different redshifts.The 
dashed dark blue line is the lognormal model. The smoothing angle 
for this set of data was 52.5 r/ . 



Fig. 5.— The theoretical model (red line) shown against the 
probability densities (blue circles) obtained from the quoted TPM 
simulation (Bode Sz Ostriker 2003) for different redshifts. The 
dashed dark blue line is the lognormal model. The smoothing 
angle for this set of data was 90" or 1.5'. 


the distribution, the lognormal does moderately well as 
long as probability densities are above ~ 10 -2 , below 
which it departs from the data points. Our model, how¬ 
ever, appears to describe these heavy tails acceptably. 
Thus the lognormal proves to be a poor approximation 
for the very high surface density tails, except for very low 
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Fig. 6. — Comparison of the lognormal and the Gaussian prob¬ 
ability density functions with the model proposed in the paper. 
The data used for this illustration is for a redshift of 0.7 and a 
smoothing angle of 52.5 // . 


or very high redshifts. In the latter limit, both of these 
models tend to the Gaussian, as will be shown below, 
and become almost indistinguishable. Another aspect of 
the distribution that the lognormal fails to describe com¬ 
pletely is its skewness. The simulated distributions fall 
sharply toward the low surface density end, whereas the 
lognormal hugely overestimates the probability density 
in this region. The proposed model appears to capture 
this feature of the distribution with appreciable accuracy. 
It goes without saying that the Gaussian fails in all re¬ 
spects, except at high redshifts and of course, for large 
enough smoothing angles for which the underlying Gaus¬ 
sian initial spectrum of perturbations remains dominant, 
(cf. Fig. 6). 

Now, we shall discuss some limiting behavior of the 
above model. For compactness, we shall drop the sub¬ 
script 2 in cr 2 in the following discussion. Let us start by 
considering the exponent in the above formula, namely, 


(In a; + w 2 /2) 2 (l + A/x) 
2uA 


(15) 


Since, u> 2 and A appear to be well behaved functions of 
cr 2 , one can write, 


UJ 2 = CoCT 2 + 0(cr 4 ) 

A = ao + aid 2 + 0(a 4 ) (16) 


We can expand each component of equation (15) as, 

(lnx + u> 2 /2) 2 = (lnx) 2 + co(lnx)cr 2 + 0(a 4 ) (17) 
1 A ~ = 1 T ——■ T —-cr 2 + 0(d 4 ) (18) 

X XX 

Therefore, the argument of the exponential can be ex¬ 
panded as, 

(lnx + w 2 /2) 2 (1 + A/x) 


1 


2 w 2 


2cocr 2 

(lnx) 2 


(lnx) 2 ( l + —) 
x 


a i- 


a-o , 


+ 0(a 4 ) 


+ Co(l + —) Inx j cr 2 

(19) 


For cr 2 <C 1 up to order 0(cr 2 ), this reduces to, 


(lnx + w 2 /2) 2 (1 + A/x) 

2U 2 


1 

2cqCt 2 


(lnx) 2 (l + — )(20) 
x 


For small cr 2 , the range of x over which the exponent 
remains appreciable, is also small, 

x =1 + 8 


where S 2 ~ 0(cr 2 ). In this approximation, 
lnx = ln(l + S) ~ J 


and 


- = 1-5. 

x 


5 2 


This implies, 

(lnx + w 2 /2) 2 (1 + A/x) 

2 to 2 2co<7 2 

Under these approximations the model reduces to, 


-(1 + ao) (21) 


N 

— exp 


(lnx + w 2 /2) 2 (1 + A/x) 

2J 2 


~ N exp 
= IV exp 


5 2 


2co cr 2 

5 2 

2/3 2 cr 2 


(1 


Where we have defined 0 = (1 + ao)/co- One can 
easily verify that the condition in equation (14) gives, 
cr 2 = cr 2 a 2 , implying a = 1. So the normalization condi¬ 
tion, equation (12) gives us N = l/(\02ncr). This is also 
graphically illustrated in Fig. 7. Thus our model cor¬ 
rectly reduces to a Gaussian for small values of variance. 

For small cr, we have 1 + a D = Co — 3 and therefore, 

< jJ 2 ~ 3cr 2 + 0(cr 4 ) 


At this point, it is important to note that the surface 



Fig. 7. — At small a 2 , the parameters (1 + A) and become 
identical implying (1 + ao) = Co — 3 (cf. equation (16)). Also, the 
normalization parameter N tends to the Gaussian value 1 / . 

These imply that the distribution tends to the Gaussian form with 
variance a 2 for small values of the latter. 

density, E, that we have been considering so far, con¬ 
tains in it the contribution from the smooth background 
matter density (E). However, lensing is sensitive only 
to the departures in the surface mass density from the 
average rather than its total value (Petters et al. 2001, 
equations 3.32 and 3.48). One should note in this con¬ 
nection that the values of E relevant to strong lensing 
events are usually much higher than the background den¬ 
sity (E). This often leads to neglecting the subtraction of 
background density in the definition of k when discussing 
strong lensing. In many cases, this is indeed a very small 
effect. However, we have found appreciable differences, 
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especially in the multiple plane lensing probabilities (see 
Section 4), if this correction is not taken into account. 
Therefore, we define the lensing surface mass density £;, 
that usually appears in the lensing equations as, 

£, = £-<£>. (23) 


For the purpose of studying the possibility of strong lens¬ 
ing by a plane, it is important to notice that a multi¬ 
ple image of the source results if the beam encounters a 
background-subtracted surface density, £; along its path 
which is in excess of the critical surface mass density 
^crit^) at that redshift, which is given by, 


£crit (z) 


c 2 D s 
4 t:G A, A ’ 


(24) 


lognormal form studied by Taruya et al. (2002) in de¬ 
scribing the line-of-sight or effective convergence, given 
by (see, e.g. Bartelmann & Schneider 2001), 

J dxg(x,Xs)8{x), (33) 

where \ is the comoving distance, \s is the comoving dis¬ 
tance to the source, <5 = Ap/ ( p ) is the density contrast, 
and 

9(x,Xs) =r(x)— 7 — 7^(1 + z), (34) 

r{Xs) 

r(x) being the comoving angular diameter distance. This 
will be the subject of a future paper. 

4. STRONG LENSING PROBABILITIES 


where D s , Di and Di s are the angular diameter distances 
from the observer to the source, the observer to the lens, 
and the lens to the source, respectively. 

In order to address the problem in which one, two or 
more planes together become supercritical, it is conve¬ 
nient to change the variable of interest from x = to 

the convergence, 

k =^- = ^L{x- 1). (25) 

^-'crit ^crit 

The PDF g(n) for the convergence is related to that of 
x, f(x) by a simple Jacobian transformation, 

g(n) = af(an + 1), (26) 

with, a = 

A physical meaning to the variable a can be attached 
by defining the minimum value K m in of the convergence 
k 1 corresponding to an “empty” beam, with £ = 0. In 
this case, equation (23) implies £/ = — (£) and hence, 


A 'min — 

^crit OL 


(27) 


Therefore, one may write a = —-—, and the PDF for the 
convergence can be rewritten as, 


ff(«0 



Kmin 



or, written out in its entirety as, 


(28) 


g(n)dK = N 

x exp 

x 


(ln(l + Ac/| Kmin \)+u 2 /2) 2 {1+A/{1 +k/\ ^min ill 
2^- 


dhz 


^“l - 1 ^min | 


(29) 


The constraint equations, (12), (13) and (14), translate 
to, 

Normalization : g(n)dn = 1, (30) 

By defn. (k) = 0 : Kg(n)dK = 0, (31) 

Observed variance: n 2 g(k)dn = (k 2 ) . (32) 


This transformation leads us to the possibility of expand¬ 
ing the applicability of our model beyond the thin mass- 
sheet case that we have been considering so far. Although 
we will not consider this here, it may be interesting to 
ask how much improvement this model provides over the 


The necessary condition for a strong lensing situation 
is K e ~ K i A 1, where i is the index of a plane in 
front of the source. As the convergence k at a point on 
each plane can be either positive or negative, a super¬ 
critical plane does not necessarily imply a supercritical 
line of sight . In fact, as Premadi & Martel (2004) note, 
underdense regions tend to occur around cluster lenses. 
However, when a single plane or a number of planes, col¬ 
lectively, is supercritical, the probability that negative 
convergences on some of the remaining planes will make 
the line of sight subcritical is actually small. Neglecting 
this real but small correction for the moment, we shall 
show that there exists a possible analytical calculation 
to compute the probability that a single plane produces 
strong lensing. Then, in a following subsection we shall 
turn to a more accurate numerical (Monte Carlo) ap¬ 
proach. 

Notation: Let T denote the probability that the i-th 
plane is supercritical, i.e. 

/ OO 

g(K(zi))dK(zi). 

(35) 

4.1. Probability that one and only one plane is 
supercritical 

This is equal to the probability that the first plane 
is supercritical and all the other planes are not, or the 
second plane is supercritical and the other planes are not 
, and so on: 

N N 

'P one = d~i (1 ~ Fj) 

*=1 j = 1 

3 i 
N 

( 36 ) 

*= i 

where N is the number of planes in front of the source, 
and the last step follows from the fact that Ti <C 1, Vi. 
An extension of this scheme for two or more planes be¬ 
comes increasingly cumbersome and impractical as the 
number of planes increases. Also, as we already pointed 
out, this method does not correct for the fact that some 
lines of sight may become subcritical due to negative con¬ 
vergences on the planes other than those considered in 
the probability calculation. Therefore, we turn to the 
more efficient Monte Carlo approach to study these prob¬ 
abilities. 
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Fig. 8. — Top: The probability, as a function of source redshift, 
that a single plane alone renders the line of sight supercritical. 
The smoothing angle for this figure was 7.5 // . The solid line is 
the analytical approximation (cf. equation (36)). The filled cir¬ 
cles represent the single-plane probabilities generated by a Monte 
Carlo simulation which does not take into account the correction 
discussed in the text. The open circles represent the same, but 
with the correction applied. As discussed in the text, the analyt¬ 
ical calculation slightly overestimates the true probabilities (open 
circles) as it includes rare cases for which, although the single plane 
is supercritical, the line of sight is rendered subcritical by negative 
convergence on some other plane or planes. The open circles are 
joined by a spline (dashed line) as a guide to the eye. Bottom: 
The probabilities, shown as the percentage of the corresponding 
single-plane values, that two, three, or four planes together render 
the line of sight supercritical. The points are from the Monte Carlo 
study and are joined by splines (dashed lines) to aide the eye. 


4.2. Monte Carlo method. 

We randomly draw values, /«,, from the PDFs of k on 
the planes lying in front of the source and consider only 
those cases for which the total convergence is greater 
than unity. With k,; > 1 satisfied, if the highest value 
among the Kj’s is greater than unity we identify this as 
a single plane strong lensing case. If this is fails, we 
consider the two highest values of k from our crop and 
look at the sum of the two. If the sum exceeds unity then 
we identify this as a two plane lensing case. We extend 
this method to include three or more planes. 

4.3. Results 

As we have already shown, the probabilities for the 
single-plane cases are easy to compute analytically (cf. 
equation (36)), as long as we do not care whether un¬ 
derdensities on the other planes are making the line of 
sight subcritical. With this assumption, the analytical 
calculation slightly overestimates the probabilities as ex¬ 
pected. However, it is useful as a quick and handy way 
to get an upper limit to the probabilities. The results 
from the Monte Carlo study do take the effect of under- 
dense planes into account. The results from both the 
analytical and the Monte Carlo approach, for a smooth¬ 


ing angle of 7.5", are displayed in Fig 8. We note from 
the figure that when the constraint y) k > 1 is not ap¬ 
plied in the Monte Carlo study and cases with « max > 1 
but K < 1 are wrongly counted as legitimate single¬ 
plane cases, the single-plane probabilities obtained ana¬ 
lytically match those from the Monte Carlo simulation, 
as expected. When the correction is taken into account, 
the points from the Monte Carlo study lie slightly below 
the analytical values, because we lose a few cases due to 
the fact that the line of sight was subcritical although the 
single plane was supercritical. We notice that for source 
redsliifts, z s ~ 1.0, the contribution of secondary clumps 
of matter to strong lensing is almost negligible (~ 2 — 3% 
of the single-plane cases). However, for larger redshifts 
secondary planes become more and more important and 
contribute up to about 15% of the single-plane cases for 
z s = 6.0. Our results are somewhat lower than those ob¬ 
tained in a previous study by Wambsganss et al. (2004b). 
This may be due to various reasons, including the fact 
that we have used a smoothed surface density distribu¬ 
tion and are considering a larger opening angle than was 
done in Wambsganss et al. (2004b) . However, our results 
seem to be in good agreement with those quoted in Dalai 
et al. (2005), who report a secondary plane contribution 
of < 7% for a source redshift, = 2.0 and angular scales 
> 10". 

5. ERROR IN THE MASS ESTIMATION OF THE MAIN 
LENS 


As discussed in the previous section, there is a finite 
and non-negligible probability that more than one lens 
is responsible for causing a strong lensing event. We 
may ask, by how much we can underestimate the mass 
in lensing cluster or galaxy by assuming that it is the 
only object responsible for lensing. We study this by 
again shooting rays through the planes and considering 
cases where the total convergence along the line of sight, 
Ktot = K i greater than unity. Within these cases, 
we look at the difference between K to t an d the maximum 
value, K max among the k,:’s on the planes in front of the 
source. This maximum value would be the true value of 
the convergence due to the primary cluster which is the 
major contributor to the lensing. If we think that the 
other planes do not contribute significantly, we would 
assign the value Ktot to it instead of K max and this will 
be an error because statistically the contribution due to 
secondary planes is not negligible. To obtain an estimate 
of the error, we repeat the exercise several times and 
draw up the probability distribution of the percentage 
error in k. defined as, 


^tot ^max nr. 

e =- x 100 


(37) 


^max 

We do this for a source redshifts of 2.0 and 3.0 and for 
each redshift, draw up two probability distributions of e 
corresponding to the smoothing scales of 7.5" and 30". 
The results are displayed in Fig. 9. The distributions for 
the smaller smoothing length (~ 7.5") has a much longer 
tail (not shown), extending to very high values of e, than 
that for the larger smoothing angle. This is expected 
because the probability density function of k has a heav¬ 
ier tail for small smoothing angles, than it has for larger 
ones. We have also attempted to quantify the central 
value and the dispersion in e by fitting a Lorentzian-like 
function to the data, given by, 
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TABLE 1 

Parameters obtained by fitting the formula in equation 

(38) TO THE PDF OF e. 


Source redshift (z s ) 

Peak centroid (ai) 

HWHM (a 2 ) 


e = 5" 0 = 20" 

6 = 5" 

II 

to 

o 

2.0 

1.52 4.82 

4.93 

4.39 

3.0 

1.45 4.55 

6.52 

4.62 


/ 0 ) 


_Oo_ 

\{x - ai)/a 2 ] 2 + 1 3 


(38) 


where the best fit values of the parameters a\ and a 2 
are interpreted as the centroid and the dispersion (half¬ 
width at half-maximum, or HWHM) of the distribution. 
These values are tabulated in Table 1. 

From the table we note that the dispersions in the val¬ 
ues of e around the peak increases with increasing source 
redslrift for a fixed smoothing angle. This is expected 
because for higher redshifts there are more intervening 
planes and hence a larger range of masses in clumps lin¬ 
ing up with the main lens. Also, at a given redshift, the 
spread in the values of e is larger for a smaller smoothing 
angle. This, basically reflects the fact that the variance 
of k is higher for smaller smoothing angles. We also note, 
that although the distributions for the smaller smooth¬ 
ing angle has a much longer tail (not shown), the most 
probable error corresponding to it is lower than that for 
the larger smoothing angle. This basically means that 
although multiple-plane cases are less probable for the 
larger smoothing angle, when they happen, they carry 
a larger contribution of mass from the auxiliary planes, 
than in a typical small angle case. For instance, in a two- 
plane case the value of k on the secondary plane typically 
come from the region around the peak of the PDF of n 
on that plane. For a small smoothing angle, the peak is 
close to k < 0 (the ray mostly hits slightly underdense 
regions), whereas as the smoothing angle grows the peak 
moves towards more positive values of n. So in a su¬ 
percritical line of sight, secondary planes will typically 
contribute higher values of k when the smoothing angle 
is larger. This also applies to all the planes along the line 
of sight in general, so that for a larger smoothing angle 
the residual sum from all other planes, K t ot — Kmax, will 
be higher than that for a smaller smoothing angle. 


6. EFFECT OF INHOMOGENEITIES ON ANGULAR 
DIAMETER DISTANCES 

Inhomogeneities in the distribution of matter induce 
fluctuations in the measured value of the angular diam¬ 
eter distance at a fixed redshift, z. The problem has 
been attacked directly from general relativity (Zel’Dovich 
1964; Refsdal 1970; Dyer & Roeder 1974; Sasaki 1987; Fu- 
tamase & Sasaki 1989; Watanabe & Tomita 1990; Segal 
1993; Kantowski et al. 1995; Linder 1998; Kantowski et 
al. 2000; Mortsell 2002; Demianski et al. 2003) producing 
results of varied complexity. 

In the following, we shall show that under certain sim¬ 
plifying assumptions, the differential equation for an¬ 
gular diameter distances can be obtained as a contin¬ 
uum limit to multiple gravitational lensing by successive 



Fig. 9. — Probability distribution of the percentage error, e, in 
the mass of the main lens for two different smoothing angles. Top: 
Source redshift, z s = 2.0. Bottom-. Source redshift, z s = 3.0. 


mass-sheets. This will provide a framework for study¬ 
ing angular-diameter distances in an inhomogeneous uni¬ 
verse as a multiple deflection problem. The calculations 
are much in the spirit of Hadrovic & Binney (1997), al¬ 
though we have generalized their results to include a non¬ 
zero cosmological constant. 

Let us use the convention that D and D 1 respectively, 
denote angular-diameter distance before (empty beam) 
and after lensing is taken into account. 

Fig. 10 shows a light ray that passes two mass sheets at 
impact parameters (measured with respect to the fiducial 
ray, OO') and £ 2 an d is deflected through angles a.\ 
and ct 2 , respectively. A simple geometrical argument 
(see, for e.g. Hadrovic & Binney 1997) shows that the 
impact parameter, on the j-th plane is related to ^ 
through the relation, 

D _ 

= 1 ~~ El Dij OL i(£,i)i (39) 

U\ , 

1 ?.= I 


where Di is the enrpty-beam angular diameter distance 
to the i-th plane and Dij denotes the enrpty-beam an¬ 
gular diameter distance between the i-th and the j -th 
planes. 

If we assume that the mass distribution inside the 
beam on each mass sheet is uniform, with a surface den¬ 
sity £, then the angle of deflection is given by, 


a = 



(40) 


Hence in this case the impact parameters are all parallel 
and satisfy 





4t tG 


j -1 




(41) 
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Fig. 10.— Lensing by a series of mass sheets. 


The angular diameter distance to an object with diame- 
ter is = £*/(£i/Di). Therefore, on taking the 

scalar form of equation (41) and dividing throughout by 

(£i/Di) we find that the equation satisfied by the true 
angular diameter distances Dj has the form, 

_ 4 g i _ 

Dj = Dj — ^2 ^ ' Ej D, D, :/ . (42) 


6.1. Equation for D{y,z). 


Let us denote the redshifts of the i-th and the j-th 
planes as y and 2 , respectively. Also, let the continuous 
version of the distance Dy be denoted by D(y, z ). 

D(y, z) is the angular diameter distance in an empty 
beam (vanishing Ricci tensor and shear) and hence the 
focusing theorem (Schneider et al. 1992, equation (3.64)) 
takes the form, 


d 2 D 

d^ 


= 0, 


(43) 


where r is the affine parameter, defined through, 


dr = 


_ dz _ 

(1 + z) 2 fl m ( 1 + z)^ + D&(1 + z) 2 + Da 


where we have defined, 


dz 

wr 

(44) 


g(z) = (1 + z) 2 \JQ .„, (1 + z) 3 + D&( 1 + z) 2 + Da- 

With the help of equation (44), equation (43) can be 
expressed in terms of redshift as, 

9( z )§: = 0. (45) 


The initial condition is given by (Schneider et al. 1992, 
equation (4.53), same arguments, but with Da ^ 0), 


dD (y,z) = (1 + y) 

dz z=y y(y) 


(46) 


6.2. Transition to continuum 

Now we shall fill the beam with a cosmological matter 
density and use equation (42) to calculate the angular- 
diameter distance to a redshift z. We consider matter 
in the beam to be sliced up into discs of infinitesimal 
thickness. The disc that lies between redshifts y + dy 
and y has surface density 

dE = p m (y)%Ldy, (47 ) 

dy 


where, the proper distance, ?’ pr0 p is given by, 

dy 


dr prop — cdt 


Ho (1 + y)\/D m ( 1 + y) 3 + Dfc( 1 
c (1 + y) 


V? 


Da 

(48) 


Ho g(y) 

and the physical density at redshift y, p m (y) is given by, 

Pm(y) = ^n m [l + 6(y)](l + y) 3 , (49) 

where D m is the dimensionless density parameter in the 
present universe and we have entertained the possibility 
of having inhomogeneous distribution of matter through 
the perturbation J(y) in density. 

Using equations (48), (49) and (47), we have, 

dE = p m [i + <5(y)](i + y) 3 4r^-z7~p- d y 


Ho g{y) 


3cH 0 
87 tG 


{1 + 6(y)} 


(! + y) 
g(y) 


41 


dy. (50) 


The continuum limit of the multiple deflection equation 
(42), after making the replacements: Dj — > D(z) 1 Di —> 
D(y ), Dij —> D(y, z) and J2i=i —> Jo dy, takes tke f° rm , 


D(z) = D{z)~ 


J dy ~^J D ( y ^ D ( y ’ z J 


(51) 

Using equation (50) and expressing all distances in units 
of c/Ho , the above equation becomes, 


D(z) = D(z) 


— s d 


dy 


I 1 + <%)](!+ y) 4 


D{y)D{y,z\ 52) 


io g(y) 

We can convert this integral equation to a differential 
one as follows. Upon differentiating equation (52) with 
respect to z, we have, 

dD _ dD 
dz dz 

,(1 


— - 

2 


dy [! + %)]■ 




and writing D m [l + <5(y)] as D m (y) for brevity, we further 
have, 


(; k S(y ’ z) 



dz J 

-|D m (z)(l + z) 4 D(x) 


y=Z 


— 3 Q 
2 ^ L r\ 


dy 


(1 -)- yU 

{i + %)} L 44 L D(y) 


g(y) 


d 


dz 


d 


x AT \a{z)^-D{y,z) 


dz 


The first and the last terms on the right hand side of the 
above equation vanishes by virtue of (45). Using (46), 
we finally have the differential equations for D{z) as, 

9( z )jZ + I ^(f + S(z)) (1 + zfD(z) = 0. 

(53) 


dz 
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Fig. 11.— Top: The mean angular diameter distances obtained 
by Monte Carlo sampling of surface densities from the lens planes. 
The continuous line represents the smooth Friedman Robertson 
Walker (FRW) angular diameter distance as a function of redshift. 
The data are for an angular smoothing scale of 7.5". Bottom: 
The ratio, cr^(z)/(D(z)), of the standard deviation to the mean 
of the angular diameter distance to each redshift, for two angular 
smoothing scales, 7.5" and 30.0". 


If <5 = 0 then equation (53) is the correct differential 
equation for angular diameter distance in an FRW uni¬ 
verse. If the ray propagates through a smooth under- 
dense region having a constant comoving density 
(0 < a < 1), then the equation reduces to the historical 
Dyer-Roeder equation. 

Hence, we have shown that the angular diameter dis¬ 
tances can be studied as multiple deflection lensing prob¬ 
lem. It is important to discuss the simplifying assump¬ 
tions that lead to this derivation. Firstly, we have as¬ 
sumed that the focusing of rays propagating over cos¬ 
mological distances is unaffected by shear due to clumps 
of matter lying outside the beam. Futamase & Sasaki 
(1989) and Watanabe & Sasaki (1990) show that shear 
contributes negligibly to focusing as long as the scale of 
inhomogeneities are greater than or comparable to galac¬ 
tic scales. Since we study this problem with smoothing 
scales ~ 10", this assumptions seems to be a reasonable 
one. Secondly, we have assumed that matter that lies 
between the redshifts z and z + dz forms a uniform disc. 
This is a nontrivial approximation but one that is par¬ 
ticularly suited to our setup, in which we are looking 
at surface mass distributions on planes which have been 
already smoothed by a window of a certain angular size. 

6.3. Implementation and results 

For each beam, we approximate the continuous distri¬ 
bution with a Monte Carlo approach utilizing 38 values 
of S g from the PDFs corresponding to the 38 planes and 
compute the angular diameter distances out to various 
redshifts using equation (42). In the top panel Fig. 11, 
we show the result of a Monte Carlo run with 3 x 10 4 


Fig. 12.— The probability distribution function of the angular 
diameter distance to four different redshifts, z = 1.34, 2.01, 2.99 
and 4.49 and for two smoothing angles 6 = 7.5" (blue histograms) 
and 30" (red histograms). The dashed vertical lines mark the 
position of the average or FRW angular diameter distance at each 
redshift. 


beams, each having an opening angle of 7.5". As ex¬ 
pected, the mean angular diameter distances perfectly 
match the smooth Friedman Robertson Walker (FRW) 
values. However, as displayed in the bottom panel of the 
same figure, the ratio of the standard deviation to the 
mean value of the angular diameter distance, increases 
with increasing redshift, amounting to a probable error 
of ~ 5% at a redshift of ~ 2.0 and rising up to ~ 9% for 
a redshift ~ 7.0, for the smoothing angle of 7.5" . The 
dispersions diminish with an increase in the smoothing 
angle, as is evident from the data points for the smooth¬ 
ing angle of 30". This is expected because with a larger 
smoothing angle, the surface densities encountered by the 
beam departs from the mean density with lesser proba¬ 
bility. In Fig 12, we show the probability distributions of 
the angular diameter distance D, at four redshifts. These 
distributions are skewed towards values of D lower than 
the FRW values and the skewness is higher for the smaller 
smoothing angle. The skewness implies that a typical 
ray more frequency passes through an underdense region 
than it does an overdense one, and as discussed before 
in relation to the error in lens mass estimation, this bias 
becomes less important as the smoothing angle increases. 
Since various combinations of the angular diameter dis¬ 
tance enter the calculations of lensing, the dispersion in 
its value, as displayed here, may systematically affect the 
use of lensing to determine of the properties of lenses or 
aspects of the underlying cosmology. 

7. CONCLUSIONS 

In the present paper, we have studied gravitational 
lensing with a multiple lens plane approach and have pro¬ 
posed a new analytic model for the dark matter conver¬ 
gence on the lens planes. Like the lognormal, this model 
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is entirely defined by (ft 2 )., the variance of ft on a given 
lens plane and for a given smoothing angle 9. We showed 
that these variances can be theoretically predicted from 
the knowledge of the nonlinear power spectrum under a 
given cosmological model. Our tests indicate that the 
proposed model is far better than the lognormal in the 
ultra-nonlinear regime of small smoothing angles (< l 7 ) 
and describes the low-probability, high k tail of the PDF 
of k with quite high accuracy. We illustrated a number 
of possible applications of the model. With the PDF of 
kappa on each plane accurately defined by our model, 
we could mimic an expensive ray shooting simulation to 
find out the importance of multiple clumps of matter in 
producing strong lensing, by a hugely faster Monte Carlo 
technique. We showed that the probabilities that a single 
plane on the line-of-sight is supercritical can be approx¬ 
imately obtained by a simple analytical calculation. We 
further demonstrated that secondary matter clumps are 
unimportant at source redshifts ~ 1 but may contribute 
up to about 15% of the single-plane cases at higher red- 
shifts. The chance accumulation of significant amounts 
of secondary matter along the line-of-sight will therefore 
result in systematic inaccuracies in the determination of 
cluster masses, typically overestimating them by about 
5% with in addition a ~ 5 — 6% standard deviation. We 
took a simplified model of the matter distribution inside 
the beam and were able to cast the problem of the deter¬ 
mination of angular diameter distances into a multiple 
lensing problem. We could then use our analytic model 
to estimate the dispersions in angular diameter distances 
due to inhomogeneities in the matter field. As a cau¬ 
tionary note, we would like to add that the for lensing 
involving small angular scales, the role of baryons be¬ 
come important and the predictions from our model in 
this regime, should therefore be looked at as some kind of 
a lower limit. An important future improvement to our 
model would therefore be the inclusion of the effect of 
baryons. It should be noted that one important advan¬ 
tage of our method over the halo model approach is the 
ability to compute strong lensing by multiple lens planes, 
which is rather difficult in the halo approach. Therefore, 
our approach provides a complementary analytic tool to 
the halo model based calculations, in the sense that our 
model will be useful to examine the effect of secondary 
matter when combined with the halo model approach. 

There are several ways in which the power of this 
model could be harnessed. As noted briefly at the end 
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of Section 3, this model may pave the way for a more 
accurate description of the line-of-sight-integrated con¬ 
vergence in the small angle regime where the lognormal 
approximation is sometimes inadequate. This issue is 
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